﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;

namespace ImageProcessing
{
	public static class Morphology
	{
		//f = image (adImage)
		//s = filter (adB)
		internal static double[, ,] Dilate(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			//f dilate s
			int nWindowOffset = (nBWidth - 1) / 2;
			int nWinTotal = nBWidth * nBHeight;
			int half = nWinTotal / 2 + 1;
			int nWinI, nWinJ;
			//double[,,] adWindow = new double[nBWidth,nBHeight,3];
			double[, ,] adOut = new double[nImageWidth, nImageHeight, 3];
			double[] adMax = new double[] { double.MinValue, double.MinValue, double.MinValue };

			unsafe
			{
				

				//calculate values
				for (int j = 0; j < nImageHeight; j++)
				{
					for (int i = 0; i < nImageWidth; i++)
					{
						if (j < nWindowOffset || j >= nImageHeight - nWindowOffset || i < nWindowOffset || i >= nImageWidth - nWindowOffset)
						{
							adMax[0] = adImage[i, j, 0];
							adMax[1] = adImage[i, j, 1];
							adMax[2] = adImage[i, j, 2];
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							//reset values
							adMax[0] = double.MinValue;
							adMax[1] = double.MinValue;
							adMax[2] = double.MinValue;
							//find max values for each channel
							for (int n = 0; n < nBHeight; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nBWidth; m++)
								{
									if (adB[m,n] >= 0)
									{
										if (adImage[nWinI, nWinJ, 0] + adB[m, n] > adMax[0])
											adMax[0] = adImage[nWinI, nWinJ, 0] + adB[m, n];
										if (adImage[nWinI, nWinJ, 1] + adB[m, n] > adMax[1])
											adMax[1] = adImage[nWinI, nWinJ, 1] + adB[m, n];
										if (adImage[nWinI, nWinJ, 2] + adB[m, n] > adMax[2])
											adMax[2] = adImage[nWinI, nWinJ, 2] + adB[m, n];	
									} //end if
									nWinI++;
								} //end for
								nWinJ++;
							} //end for
						} //end else

						adOut[i, j, 0] = adMax[0];
						adOut[i, j, 1] = adMax[1];
						adOut[i, j, 2] = adMax[2];


					} //end for
				} //end for
			}
			return adOut;
		}//end Dilate

		internal static double[, ,] Erode(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			//f erode s
			int nWindowOffset = (nBWidth - 1) / 2;
			int nWinTotal = nBWidth * nBHeight;
			int half = nWinTotal / 2 + 1;
			int nWinI, nWinJ;
			double[, ,] adOut = new double[nImageWidth, nImageHeight, 3];
			double[] adMin = new double[] { double.MinValue, double.MinValue, double.MinValue };

			unsafe
			{

				//calculate values
				for (int j = 0; j < nImageHeight; j++)
				{
					for (int i = 0; i < nImageWidth; i++)
					{
						if (j < nWindowOffset || j >= nImageHeight - nWindowOffset || i < nWindowOffset || i >= nImageWidth - nWindowOffset)
						{
							adMin[0] = adImage[i, j, 0];
							adMin[1] = adImage[i, j, 1];
							adMin[2] = adImage[i, j, 2];
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							//reset values
							adMin[0] = double.MaxValue;
							adMin[1] = double.MaxValue;
							adMin[2] = double.MaxValue;
							//find max values for each channel
							for (int n = 0; n < nBHeight; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nBWidth; m++)
								{

									if (adB[m, n] >= 0)
									{
										if (adImage[nWinI, nWinJ, 0] - adB[m, n] < adMin[0])
											adMin[0] = adImage[nWinI, nWinJ, 0] - adB[m, n];
										if (adImage[nWinI, nWinJ, 1] - adB[m, n] < adMin[1])
											adMin[1] = adImage[nWinI, nWinJ, 1] - adB[m, n];
										if (adImage[nWinI, nWinJ, 2] - adB[m, n] < adMin[2])
											adMin[2] = adImage[nWinI, nWinJ, 2] - adB[m, n];
									} //end if

									nWinI++;
								} //end for
								nWinJ++;
							} //end for



						} //end else

						adOut[i, j, 0] = adMin[0];
						adOut[i, j, 1] = adMin[1];
						adOut[i, j, 2] = adMin[2];

					} //end for
				} //end for
			}
			return adOut;
		}//end Erode
		internal static double[, ,] Open(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			//(f erode s) dilate s
			double[, ,] adOut;
			adOut = Erode(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adOut = Dilate(adOut, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			return adOut;
		}//end Open
		internal static double[, ,] Close(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			//(f dilate s) erode s
			double[, ,] adOut;
			adOut = Dilate(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adOut = Erode(adOut, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			return adOut;
		}//end Open

		internal static double[, ,] InsideContour(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			//f - (f erode s)
			double[, ,] adOut;
			adOut = Erode(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adOut =  ImageSubtract(adImage, adOut, nImageWidth, nImageHeight);
			return adOut;
		}//end InsideContour

		internal static double[, ,] OutsideContour(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			double[, ,] adOut;
			adOut = Dilate(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adOut = ImageSubtract(adOut, adImage, nImageWidth, nImageHeight);
			return adOut;
		}//end OutsideContour

		internal static double[, ,] MiddleContour(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			double[, ,] adDilated;
			double[, ,] adEroded;
			double[, ,] adOut;
			adDilated = Dilate(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adEroded = Erode(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adOut = ImageSubtract(adDilated, adEroded, nImageWidth, nImageHeight);
			return adOut;

		}//end MiddleContour

		internal static double[, ,] TopHat1(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			double[, ,] adOut;
			double[,,] adOpened;
			adOpened = Open(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adOut = ImageSubtract(adImage, adOpened, nImageWidth, nImageHeight);
			return adOut;
		}//end TopHat1

		internal static double[, ,] TopHat2(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight)
		{
			double[, ,] adOut;
			double[, ,] adClosed;
			adClosed = Close(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adOut = ImageSubtract(adClosed, adImage, nImageWidth, nImageHeight);
			return adOut;
		}//end TopHat2

		internal static double[, ,] Gradient(double[, ,] adImage, double[,] adB, int nImageWidth, int nImageHeight, int nBWidth, int nBHeight, int nNumOps)
		{
			double[, ,] adDilated;
			double[, ,] adEroded;
			double[, ,] adOut;
			adDilated = Dilate(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			adEroded = Erode(adImage, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			for (int i = 0; i < nNumOps-1; i++)
			{
				adDilated = Dilate(adDilated, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
				adEroded = Erode(adEroded, adB, nImageWidth, nImageHeight, nBWidth, nBHeight);
			} //end for
			adOut = ImageSubtract(adDilated, adEroded, nImageWidth, nImageHeight);
			return adOut;

		}//end Gradient

		private static double[, ,] ImageSubtract(double[, ,] adMinuend, double[, ,] adSubtrahend, int nImageWidth, int nImageHeight)
		{
			double[, ,] adOut = new double[nImageWidth, nImageHeight, 3];
			for (int j = 0; j < nImageHeight; j++)
			{
				for (int i = 0; i < nImageWidth; i++)
				{
					adOut[i, j, 0] = adMinuend[i, j, 0] - adSubtrahend[i, j, 0];
					adOut[i, j, 1] = adMinuend[i, j, 1] - adSubtrahend[i, j, 1];
					adOut[i, j, 2] = adMinuend[i, j, 2] - adSubtrahend[i, j, 2];
				} //end for
			} //end for

			return adOut;
		}

	}
}
